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Department of Electrical Engineering 
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Abstract 

This paper exploits the relationships among the time-domain and frequency-domain 
system norms to derive information useful for modeling and control design, given only the 
system step response data. A discussion of system and signal norms is included. The proposed 
procedures involve only simple numerical operations, such as the discrete approximation of 
derivatives and integrals, and the calculation of matrix singular values. The resulting frequency- 
domain and Hankel-operator norm approximations may be used to evaluate the accuracy of a 
given model, and to determine model corrections to decrease the modeling errors. 

1. Introduction 

Historically, system modeling has been something of an art, requiring either special 
knowledge of the system being considered, or a certain “intuition” into the modeling process 
itself. At one end of the modeling spectrum is the application of the conservation laws that 
describe many physical processes and environments. The resulting system models may be 
described by nonlinear partial differential equations in multiple spatial dimensions. These 
models can be simplified and converted into ordinary differential equations by finite-difference 
or finite-element techniques, but the resulting systems are still nonlinear and may be extremely 
complex. For example, in the field of Computational Fluid Dynamics (CFD), flow systems are 
sometimes described by literally millions of dynamic state equations. At the other end of the 
modeling spectrum, a control engineer may simply “guess” (or optimize) the parameters of a 
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low-order linear model to match the observed dynamic response of a given system. Between 
these two extreme approaches lie many techniques for creating a useful system model (Chicatelli 
and Hartley, 1997). 

Regardless of which modeling approach is used, the process always yields some 
modeling error. Still, engineers have been able to produce successful feedback control designs 
for many uncertain systems, even when deprived of accurate error descriptions or estimates. 
Apparently, the feeling that a given system model is reasonably accurate has often been 
sufficient to allow a successful control design. 

In recent years, though, the control systems community has been looking for stronger 
theoretical evidence that a feedback control designed for a nominal system model will also be 
suitable when applied to the actual system. Some such evidence has been found in the 
mathematical field of functional analysis, contributing to the evolution of the field of robust 
control (Doyle et al., 1992). The objective of robust control is to design a controller that will 
provide the desired closed-loop behavior in spite of the inevitable modeling errors and unknown 
disturbances that affect a system. To accomplish this, some knowledge of the system uncertainty 
is required. That is, the design methods that provide a priori guarantees of robust stability or 
performance invariably require prescribed bounds on the uncertainty of the system parameters or 
the system frequency response characteristics. 

The purpose of this paper is to derive useful bounds on system modeling errors from 
measured data. It will be assumed throughout that the system model is intended for control 
design or analysis; so, the derived bounds will apply to the frequency-domain descriptions, and 
be useful for the application of robust control methods derived from functional analysis. 
However, this paper is not written for the functional analys:, but for the control engineer. Many 
control engineers are trained in robust control, but few are familiar with the subtle nuances of 
functional analysis. This paper is meant to help bridge the gap from the theoretical robust 
control literature to the practicing control engineer by providing the tools necessary to quickly 
compute error bounds that will be useful in controller design. Therefore, some mathematical 
results are just stated with the appropriate references, and in-depth proofs from functional 
analysis are avoided. 
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2. Signal Norms 


A quantification of the errors in a control design model requires the measurement of the 
“size” of the error signals associated with the system. Although there are many ways to measure 
signal size (Doyle et al., 1992; Boyd and Barratt, 1991), the concept of signal norm is used here. 
The norm of a signal y(t) is generally defined as 

-i Up 

m ,= \m p * ■ (i) 

_0 

This is the general definition of the p-norm, where p is usually a positive integer. Note that y(t) 
is assumed to be defined for all t > 0. For most purposes, y(t) is assumed zero for / < 0. 

Of all the p-norms, the 1-norm, the 2-norm, and the infinity-norm are the most useful, 
mainly because of their simple intuitive interpretations. These will be discussed below. 

The 1-norm 

Setting p=l in the general definition (1) yields the definition of the 1-norm of y(f) as 

Mi (2) 

o 

Thus the 1-norm of a signal is equal to the integral of the absolute value of the signal. If, for 
example, y(t) represents the force applied over time to a body in motion, then Hy^ may represent 
the total effort or fuel expended. 

For the 1-norm of a signal to be finite, the signal must eventually decay to zero. 1 For 
signals that do not decay to zero, such as step functions or sinusoids, the integral diverges and 
the 1-norm is considered infinite. 

The 2-norm 

Setting p=2 in the general definition (1) yields the definition of the 2-norm of y(t) as 

r -|l/2 

ll>i = • (3) 

_0 

1 This statement assumes the signal is absolutely continuous, which holds for all signals produced as the outputs of 
ordinary or partial differential equations. The condition of absolute continuity is assumed throughout this paper. 
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The 2-norm is closely related to the energy of a signal. If, for example, y(t) is a voltage applied 

1 2 

is proportional to the energy delivered to the load. The 2-norm is 
also related to the variance, which is defined as 

' 1 T 

var <>) = lim -J 

L 0 

and to the RMS value (or standard deviation), which is the square-root of the variance. Note 
that, because of the division by the time interval T in (4), the variance and the RMS value may be 
finite for signals such as steps or sinusoids that do not decay. However, like the 1 -norm, the 2- 
norm is finite only for signals that decay to zero. The variance and the RMS value are 
mentioned here to provide a clearer understanding of the 2-norm, but will not be used in that 
which follows. 

The 2-norm of a signal may be defined equally well in the frequency domain as 

i 

1 “ ^2 

lY(j<»)\\ 2 = —j\Y(j(o)\ 2 d(D , (5) 

j 

where Y(j(d) is the Fourier transform of y(t), and the integration is now carried out over all 
frequencies. The factor of 2jt is for the purpose of normalization. A simple mathematical link 
between the time-domain and frequency-domain signal representations is established by 
Parseval’s theorem, which states that the signal energy in the time domain is equal to the signal 
energy in the frequency domain. This may be written in terms of the 2-norms simply as 

• (6) 

The infinity-norm 

Letting p approach infinity in the general p-norm definition (1) yields the infinity-norm 
of y(t) as 

oe 1 

ML = Jim \\y{t)\ P dt . (7) 

To 

Although this may at first appear dubious, it is resolvable. For large p, the maximum values of 
y(t) are emphasized in the integral far more than the smaller values; so, the integral is 
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approximately proportional to the maximum value of y(f) raised to the power p. Taking the 
root and letting p approach infinity yields 

||y|L=max|y ( 0|. (8) 

Thus, the infinity-norm of a signal is simply the largest value of the signal in magnitude. In 
other words, the infinity-norm of a signal represents a bound on the signal amplitude for all time. 
Unlike the 1-norm and the 2-norm, the infinity-norm of a signal may be finite even though the 
signal does not decay. 

Example 

As an example, consider the function 

y(t) = exp(-r) - exp(-2t) , t>0, 

whose graph is shown in Figure 1(a). The 1-norm of this function is calculated as 

oo j 

hi = J |exp(-0 - exp(-2/)|<* = - , 

o Z 

which is equal to the area under the graph of y(t). The 2-norm is calculated as 

i 

IM 2 = J|exp(-0-exp(-2r)| 2 d/ = ^> 

V° 

The square of the 2-norm is the area under |y(t)| 2 in the time domain, shown in Figure 1(b). 
This is equal to the area under \Y(jco)^ /2n in the frequency domain, shown in Figure 1(c). The 
infinity norm is the maximum value. 



which occurs at time t = ln(2). 

3. Two System Norms 

It is now necessary to define the norms of a system, as measures of the size of the system. 
For a system represented by a transfer function, the system norms can be defined in the 

2 To be more precise, the infinity norm should be defined using the supremum (sup), or least upper bound, rather 
than the maximum value (max) of the signal. This would take care of the case where the signal does not achieve any 
maximum value, but only approaches a certain least upper bound. 
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frequency domain similar to the signal norms in the time domain. It turns out that such a system 
norm may be interpreted as a gain of the system. Control engineers understand the gain of a 
system as a function of the input frequency, but this assumes the system inputs to be sinusoidal. 
A more general definition of the system gain is provided by the system norms (Doyle et al., 
1992; Boyd and Barratt, 1991). 

Assume that a given system is linear, time-invariant, and causal, and is described by the 
transfer function H(s). Assume also that the system is stable, so that all the poles of H(s) lie in 
the left half of the complex plane. The frequency response of the system, denoted by H(jco), is 
the transfer function H(s) evaluated on the imaginary axis, s=ja>. The system impulse response, 
denoted by h(t), is the inverse Laplace transform of H(s), or the inverse Fourier transform of 
H(jco). 

The Frequency-domain Infinity-norm 

Now the infinity-norm already defined for time-domain signals may be applied in the 
frequency domain and used as a system norm. Replacing lime with frequency in the definition 
of the infinity-norm (8) yields the new definition 

I# 0'®)L = max|//( 7 < 0 )| , (9) 

where \H(ja>)\ is the magnitude of the frequency response at the frequency ft). Thus the infinity- 
norm in the frequency domain is equal to the largest magnitude of the frequency response over 
all frequencies. Graphically, it represents the highest peak in the (Bode) magnitude plot of the 
transfer function, or the magnitude of the point on the Nyquist plot farthest from the origin in the 
complex plane. The 1-norm and the 2-norm of the frequency response are also commonly 
defined by analogy with the time-domain definitions, but these will not be considered here. 

It is clear that the frequency-domain infinity-norm of a system is the maximum possible 
system gain when considering only sinusoidal excitations. It also happens to be a useful measure 
of the system gain for inputs that are not sinusoidal. Let & general system input be denoted by 
u(t) and the system output by y(t), both of which are assumed to be scalars. (The vector case will 
be considered later.) The Fourier transforms of u(t) and y(t) are denoted as U(jco) and Y(j(o), 
respectively. With these definitions, it is well known that 

Y(j(o) = H(jco)U(j(0), (10) 
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which yields 


\Y(j(0)\ = \H(j(0)\\UU(0)\. ( 11 ) 

Observing that the frequency response magnitude |//(y'0))| at any frequency is no greater than 
the infinity-norm of the system, we obtain 

\Y (j(0)\ < \\HU<o)\\JU U g>)\ ■ ( 12 ) 

Note that the frequency-domain p-norms of Yijco) and U(jco) can be defined just as for the 
transfer function H(jco). Squaring both sides of the above inequality, integrating over frequency, 
and taking the square root yields the relation 

tyUcoi^lHUtojlpUaiy (I 3 ) 

This inequality is a bound on the output energy in the frequency domain. By Parseval’s 
Theorem, this bound also holds in the time-domain as 

U*I*</»U4. < 14) 

Thus the infinity-norm of the system frequency response provides a bound on the ratio of output 
energy to input energy in the time domain, which is a useful measure of the system gain. 

The Time-domain 1-norm 

A second measure of the system gain can be obtained directly from the time-domain 
description of the system. The time-domain equivalent of (10) is the convolution integral 

t 

y(t) = Jh(T)u(t-T)dT, (15) 

0 

where h(t) is the impulse response of the system. Taking the magnitude of both sides, we obtain 

t 

\y(t)\= jh(r)u(t-r)dr . (16) 

o 

Then the triangle inequality yields 

/ 

|y(t)| < j\h(T)u(t - t)\dt . (17) 

0 
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Since the input function u{t) is bounded by its infinity-norm (maximum value), the inequality 
(17) yields 

/ 

\y(t)\<j\h(T)\\\u\ldr. (18) 

0 

The norm of u(t) can be factored out to give 

(• \ 

\y(t)\< )\h(z)\dt |HL • (19) 

1 ° 

The integral in (19) is at its largest as t goes to infinity, at which point it approaches the 1-norm 
of the system impulse response. Consequently, the above inequality yields 

btoNMilML- ( 2 °) 

Since this bound holds for all time t, it also holds for the maximum value of y{t ) , its time- 
domain infinity-norm. Hence, (20) yields 

IML S M,IML- <2') 

Thus the 1-norm of the system impulse response is another useful measure of the system gain, as 
it bounds the maximum output value for a given maximum input value in the time domain. 

Relationship between the System Norms 

The infinity-norm of the system frequency response bounds the input-output energy ratio, 
and the 1-norm of the system impulse response bounds the input-output amplitude ratio in the 
time domain. It is now shown that, for a given system, the frequency-domain infinity-norm is 
always less than the time-domain 1-norm. 

The Fourier transform of the system impulse response is given by 

H ( ja>) = J h(t) exp (-jax)dt . (22) 

o 

Taking the magnitude of both sides and using the triangle inequality yields 

J h(t)exp(-jax)dt < J \h(t)exp(-jox)\dt = j\h(t)\\exp(-jox)\dt . (23) 

0 0 0 
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The exponential inside the integral has a magnitude of unity for all frequencies; therefore, this 
yields 

\H(jO))\<]\h(t)\dt = \\h\[. (24) 

0 

As this inequality holds for all frequencies, the left side may be replaced by the frequency- 
domain infinity-norm. This yields the final relation 

II" (HL s INI, ■ < 25 > 

Thus the frequency-domain system norm is bounded by the time-domain system norm. This 
relationship will be useful for obtaining an estimate of the system infinity-norm from time- 
domain experimental data or simulation results. 

4. Errors in Modeling 

The goal of this paper is to determine bounds on the error incurred in the modeling 

process. The modeling error can be defined simply as the difference between the approximate 

linear model and the real system or a realistic “truth model”; see Figure 2. The linear model has 
been extracted from the truth model, or from actual measured data. In general, the real system or 
truth model is nonlinear and of high order, whereas the linear model is relatively simple. For a 
given input signal, with the two systems running in parallel, an output error signal can be 
determined. 

Sources of Error 

The greatest part of the inaccuracy of the low-order linear model probably results from 
the linear approximation. Because the true system is not linear, the error depends upon the size 
of the system input. It is difficult to tell at the modeling stage how large the input can be without 
introducing unacceptably large errors. Another source of inaccuracy arises from the reduction of 
the system order, or elimination of system state variables. The determination of a low-order 
model suitable for control system design usually requires a significant reduction of order from 
the real system model, which can have millions of state variables in the case of a CFD model. 
Fortunately, the model-order reduction usually contributes relatively little to the modeling error. 
In any event, if standard model-order reduction techniques are used, the resulting errors are 
within known bounds (Glover and Partington, 1987; Glover et al., 1988). 
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Frequency-domain Representation of Error 

The system gains introduced in Section 3 are now discussed in connection with the 
modeling error and its bounds. The system model H(s) has a known frequency response H(ja>) 
that can be represented by magnitude and phase (Bode) plots or, better yet, by a polar (Nyquist) 
plot in the complex plane. At each point on the Nyquist plot, there will be some error owing to 
the inaccuracy in H(Jco) at the corresponding frequency 0). This error, denoted E(Jco), is the 
frequency response of the error system of Figure 2. The largest magnitude of the error over all 
frequencies is the frequency-domain infinity-norm of the error, represented by ||E(y'Q))|| oo . Even 
if this norm is known, it is only a bound on the magnitude of the error, and contains no 
information about the phase of the error. Consequently, all possible phase angles must be 
considered, from 0 to 2 k. Graphically, in the Nyquist plane, the error bound at each frequency is 
represented by a circle of radius ||£(y'a))||_ , centered on the Nyquist plot at that frequency. This 
means that the true frequency response of the system would lie within a “tube” in the Nyquist 
plane, as shown in Figure 3. 

This frequency-domain view of the error system presupposes that the true system is 
linear, which it is not. Nevertheless, this is the view that we shall adopt. The linear system 
norms can be thought of as providing a first-order approximation of the size or gain of a 
nonlinear system; they are useful in practice for evaluating and improving linear models of 
nonlinear systems, as the example at the end of the paper illustrates. 

Approximating the Modeling Error 

The problem now is to obtain ||E( - /©)|I . To obtain t directly would require knowing the 
error, in Figure 2, for all input sinusoidal frequencies. But it is not practical to apply sinusoidal 
inputs over a wide frequency range to either a real system or its detailed truth model. An 
alternative approach would be to make use of the results of the last section to obtain the bound 

IM. 5 HI, • <26) 

where e(t) is the response of the error system to an impulse input. If the time-domain 1-norm of 
the error system can be obtained, then it can be used as a frequency-domain error bound. 

Unfortunately, it may be difficult to obtain the impulse response of the error. This is 
particularly true if it would require an experiment with the actual physical hardware, as 
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generating and applying an impulse input may be problematical. On the other hand, accurate 
step responses are usually either known or readily obtainable from the physical hardware or 
truth-model simulations. 


Assuming again that the system is reasonably linear, the impulse response is simply the 
time-derivative of the step response multiplied by the chosen time units. This may be seen from 
Figures 4 and 5, where the operator s represents the differentiation together with the 
multiplication by the time units. The multiplication by the time units is necessary to preserve the 
overall units of the system inputs and outputs, which must be equivalent in either figure. The 
same time units must be assumed in both figures; that is, those units in which the impulse input 
is assumed to have unit area must be the same ones used to correct the system output units after 
differentiating the step response. It is safest to choose a single unit of time as the standard for the 
entire analysis, and stick to it throughout. 


A discrete-time approximation of the impulse response e{t) may be obtained from the 
step response e slep (t) via the numerical differentiation 


e(kT) = 


e slep ((k-l)T) 


(27) 


Here, the units of e(kT) are taken to be the same as those of e ste p(kT ), assuming that T is evaluated 
in the standard time units. The integral for the 1-norm of the error is then approximated by a 
sum as 

IM| 1 =Jk(o|^=X7 , K^)|. (28) 

0 *=° 

This approximation of the time-domain 1-norm of the system error provides a bound for the 
frequency-domain infinity-norm error, which can then be used on the Nyquist plot. The 
interpretation of the frequency-domain error bounds in the modeling process is now discussed. 

Referring back to Figure 3, the Nyquist plot contains the frequency response of the low- 
order linear system model, as well as the system error bound as a “tube” in the Nyquist plane. 
The radius of this “tube” represents the worst-case error magnitude over all frequencies. The 
particular frequency or frequencies for which this worst-case error occurs is not known; 
theoretically, the largest error could occur at every frequency. The low-order system model can 
be said to be relatively accurate at frequencies where its magnitude is larger than the error bound. 
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Thus, to determine a range of frequencies over which the low-order linear system model is 
relatively accurate, the error bound circles can be removed from the Nyquist plot and a single 
error bound circle drawn at the origin of the Nyquist plane, as in Figure 6. The low-order linear 
system model is then considered relatively accurate as long as its Nyquist plot stays outside the 
error bound circle at the Nyquist plane origin. At those frequencies for which the Nyquist plot 
lies inside this error bound circle, the model is considered relatively inaccurate. 

Another interesting interpretation of the error bound in the time domain comes from the 
time-domain system gain relation 

||£(0|L< \\e(t)\l |«(0|L, (29) 

where e(t) is the output error in response to an arbitrary input u(t). According to the bound, for 
inputs m( 0 that do not exceed unit amplitude, the time-domain 1-norm of the error impulse 
response is an upper bound on |e(f)| for all time. This bound assumes the error system to be 
approximately linear, and so may be most applicable for small inputs. In fact, the bound is 
usually quite conservative; so, it may still be usable for larger inputs, despite the nonlinearity of 
the real system model. 

Error Incurred in Further Reduction 

It is sometimes desirable to further simplify a model by a model-order reduction after the 
initial linear approximation. In this situation, the frequency- domain error bounds simply add, as 

total UCO)\l — ll^/frMronzanorj reduction l o-)L + IK— <HL+- < 30 > 

This means that an approximation of an approximation increases the error bound by the 
magnitude of the additional error bound. This approach stands as an alternative to recalculating 
the error bound by applying Equations (27-28) to the step response error of the newly reduced- 
order system. 

Modeling Errors in Control Design 

The field of robust control is concerned largely with accounting for unknown modeling 
errors in the control design process (Doyle et al., 1992). A simple discussion of the effect of 
modeling errors on the feedback control system design is now presented. The closed-loop 
system is assumed to be of the form given in the block diagram of Figure 7. In the diagram, H(s) 
is the low-order linear system model, E(s) is the modeling error, and the sum H(s)+E(s) 
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represents the real system. Neither the modeling error E(s) nor its infinity-norm ||i?( i /G>)|| >o is 
precisely known, but the time-domain 1-norm ||4 can be found, and is a bound on the infinity- 
norm. Therefore, the worst-case error at any frequency may be represented as 

c/®)-l4 ex pO' 0 )’ < 31 ) 

where the phase 0 of the error is unknown. The expression (31) represents the error circles of 
radius ||e|| in the Nyquist plane. 

To determine the stability of the closed-loop system, the frequency response of the loop 
transfer function is plotted in the Nyquist plane. The loop transfer function frequency response 
is given by 

T{jOi) = G(jco)[H(jco) + E(j(o)], (32) 

where G(s) is the compensator transfer function. To represent the largest possible deviation of 
T(jo)) from its nominal value, we replace E(jco) by the worst-case error, to obtain 

T{ja) = G(jco)[H(jco) + 1|4 exp(y0)]. (33) 

This expression yields the worst-case loop transfer function as 

T(jo ) ) = G(jco)H(ja>) + |GO'm)|||4 exp ijO) ■ ( 34 ) 

The last term in (34) represents the uncertainty in the loop transfer function. Note that G(jO ) ) has 
been replaced by |G(yft>)|, since the arbitrary angle 6 already accounts for the compensator 
phase. The uncertainty is represented in the Nyquist plane by error-bound circles whose radii at 
each frequency are given by |G(y‘to)|||4 . Thus, around the Nyquist plot there is a tube whose 
radius changes with frequency, as shown in Figure 8. To guarantee stability, the compensator 
must be designed so that this tube avoids the critical point — 1+yO in the Nyquist plane. 

In some cases, the analysis may be simplified somewhat by considering the frequency- 
domain infinity-norm of the controller. Replacing |G(y'<J))| by its largest value ||G(ya>)|| M yields 
the worst-case loop transfer function 

T(j(0) = G(jco)H (ja>) + ||G(y'G))||J|4 exp(y0) . (35) 

Using this expression, the Nyquist plot of the nominal compensated plant can be plotted with 
error bound circles of equal radii, as shown in Figure 9. The radius of each circle is the product 
of the compensator and error norms. If the resulting tube in the Nyquist plane misses the critical 
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point — 1+/0, then the system is guaranteed to be stable for the given uncertainty of the plant. An 
equivalent condition is shown in Figure 10, where the error bound circles have been removed 
from the frequency response, and a single circle has been drawn centered at the critical point in 
the Nyquist plane. Thus, the major effect of the uncertainty error bound is to effectively enlarge 
the critical point in the Nyquist plane to a circle. This circle must be avoided by the closed-loop 
frequency response to guarantee stability in the presence of the bounded modeling errors. 

In many control systems, such as those that require integral control for tracking, the 
compensator has a very large infinity-norm because of a pole on or near the jco axis. In this case, 
the tube in the Nyquist plane would consist of a very large circle at each point of the Nyquist 
plot, or equivalently at the critical point in the Nyquist plane. Then the simplified condition 
derived from (35), while still valid, is not particularly useful, and the circles of unequal radii 
derived from (34) should be used. 

5. Norms Based on the Hankel Operator 

This section presents a few more system norms based on the definition of the Hankel 
operator of the system (Partington, 1988). If the impulse response of a system is denoted by hit), 
the Hankel operator T of the system is defined by the relation 

(n<X0 =\h{t + x) h(-t)<2t , f >0. (36) 

o 

This is an integral operator mapping all past inputs u into future outputs Tu. That is, the output 
Tm of the Hankel operator is the same as the output of the system for t > 0 in response to an input 
u applied only for t < 0. 

The Hankel Singular Values 

The Hankel operator turns out to be rather easy to work with, and also useful in finding 
reduced-order models for complex systems. To this end, we will write the expressions that 
define the eigenvalues and eigenvectors of the Hankel operator. The n * eigenvector of the 
Hankel operator is an input function u„(t) defined for t < 0 such that 

0 n u n (-t)=]h(t+T)u n (-T)dT , t> 0, (37) 

0 
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where the scalar cr„ is the corresponding eigenvalue. The eigenvector relation can equally well 
be written as 

<7 n v n (T)=jh(t + T)v n (t)dr , T>0, (38) 

0 

where v„(r) is defined for t > 0. The functions u„(t) and v„(/) are essentially the same, except that 
one is defined for positive t and the other for negative t. 

Because the kernel h(t+f) of the Hankel operator depends on both t and r in the same 
way, the operator is called “self-adjoint,” which means that it has many of the properties of a 
symmetric matrix. Therefore, it turns out that the eigenvectors of the Hankel operator are 
equivalent to its singular vectors, and are mutually orthogonal. The eigenvector relations for the 
Hankel operator may therefore be written in the form of singular vector relations as 


(t)=^h(t + T)u n (-T)dz , t> 0, 

0 

(39) 

(-T) = ^h(t+T)v n (t)dt, T>0, 

(40) 


0 


where the functions u n and v„ are known, respectively, as the / 1 th left and right Hankel singular 
vectors, and u n (-t) = v n (t) for t > 0. The Hankel singular vectors are also referred to as the 
Schmidt pairs of the Hankel operator. The scalar a n is known as the / 1 th Hankel singular value of 
the system. The number of nonzero Hankel singular values (and corresponding Schmidt pairs) is 
equal to the order of the system. The operator describing a spatially distributed system may have 
infinitely many nonzero Hankel singular values. 

The Hankel singular values of a real system are always real and positive, and are usually 
listed in decreasing order so that Ci is the largest. They are related to the gain of the system, and 
are used to define several system norms, as follows. 

The Hankel norm of a system with the transfer function H(s), denoted \\H (s)| w , is 
defined to be the largest Hankel singular value; that is, 

< 41 > 
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This norm represents the maximum ratio of future output energy to past input energy [Boyd and 
Barratt]. In mathematical terms, one would say that the Hankel norm of the system is the 
induced 2-norm of the Hankel operator, or 

|«( J )|„=|r|| = m_axli = ff , (42) 

Another norm of H(s ) related to the Hankel operator is the nuclear Hankel norm, denoted 
\\H(s)\\ n h . It is defined as the sum of the Hankel singular values. That is, 

l««IL„ = 5>. ■ <«> 

n = 1 

The Hankel operator is said to be nuclear if its nuclear norm is finite. 

Yet another norm of H(s ) related to the Hankel operator is the Hilbert-Schmidt Hankel 
norm, denoted ||#( 5 )|| W _ 5 _ H • ft is defined as the square root of the sum of the squares of the 
Hankel singular values. That is, 

• <*» 

( n=1 J 

This can be thought of as the Frobenius norm of the Hankel operator, and is also referred to as 
the Frobenius-Hankel norm, denoted ||//(s)|| • If this norm is finite, the Hankel operator is 

said to be Hilbert-Schmidt. Interestingly enough, the square of the Hilbert-Schmidt Hankel norm 
is equal to the time-weighted impulse response energy, 

< 45 > 

n - 1 0 

It can be shown that this quantity is also equal to the irea inside the Nyquist plot of the 
frequency response of the system (Hanzon, 1992). 

The above system norms are known (Glover et al., 1988) to be related by the inequalities 

l«<4, < |wo«ol £ M, s 2i "MIL • < 46 > 
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The Hankel Approximation of a System 


If the Schmidt pairs of an integral operator are known, the kernel of the operator is given 
by an expansion similar to the singular-value decomposition of a matrix. In the case of the 
Hankel operator, the kernel is the system impulse response, and is expanded as 

fc(f + T) = jTff B v„(f )«„(-*). (47) 

n=l 

This is called the Hilbert-Schmidt expansion of the Hankel operator. Given this expansion, an 
approximation of the impulse response may be found by truncating the series as 

+ = (48) 

71=1 


where q is the order of the approximation. A system with the impulse response h(t) will 
constitute a ^ th -order model of the original system with impulse response h(t ) . Based on this 

idea, Glover et al. (1988) derived a state-variable description of this g^-order Hankel 
approximation, solely from the first q Hankel singular values and the Schmidt pairs of the 
original system. The resulting state- variable parameters are 


B i =<T i v l .( 0) 

C, =«,.(' 0) 


A, =— ~ m , 2 (°) 

B]Bj -ofC'tCj 


(49) 

(50) 

(51) 

(52) 


where the indices i and j are bounded by q, the order of the approximation. The given state- 
variable parameters describe a balanced realization in output normal form. Notice that the 
Schmidt pairs are evaluated at time t = 0. It is assumed that the Schmidt pairs are normalized, so 
that |m ( || 2 = ||v ( |j 2 = 1 for each i. 

The accuracy of the Hankel approximation is measured by the truncated Hankel singular 
values. Specifically, if H(s) describes the original system and H(s) the q ^- order Hankel 
approximation, the Hankel norm of the approximation error is given by 

|tf(.s)-H(*)|| w =O q+l , (53) 
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where o q + 1 is the largest Hankel singular value of the original system that has not been taken into 
account in the approximation. Further, a frequency-domain bound is given by 

|tf(*)-£(s)|L<2 jj/r,. (54) 

n=q + 1 


Application to the Error System 

If h(t ) represents the impulse response of the error system E(s), and the associated 
Schmidt pairs and Hankel singular values can be found, then the upper and lower bounds on the 
frequency-domain error norm are given as 

|£( J )|„<||£<y<»)|L<2|£( S )|„_ ( , . (55) 

Furthermore, the formulas (49-52) may be used to determine a system E(s) of order q that 
approximates the error system. Recall that the actual error system is not precisely known, and 
may be of extremely high order. The g^-order approximate system E(s) is used as a correction; 
that is, it may be added in parallel to the system model to decrease the error. The error 
remaining after the Hankel correction is given as E(s) = E(s) - E(s ) , and is bounded in the 
frequency domain by 

|£(-HL * 2 i>* ’ < 56 ) 

Tt=<7+1 

as well as by 

||£0'a»|L - HO ~ ? (0|| » ( 57 ) 

where e(t) is the impulse response of the correction E(s) . 

Approximation of the Hankel Operator 

So far, we have assumed that the Hankel singular values and singular vectors are known. 
However, finding them usually requires a linear system description in state-variable form. 
Lacking this, we turn to numerical methods for approximating the calculations of Glover et al. 
(1988). 

The Fredholm approach to solving integral equations (Lovitt, 1950) is adapted for 
obtaining a useful approximation. Basically, the Fredholm approach is to divide the integration 
limits for the Hankel operator into evenly spaced intervals in t and T, and then to add up the 
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resulting samples of the kernel using the rectangular rule for integration. The integral 
relationship for the eigenvalues and eigenvectors, 


c n u n (t) = jh(t + T)u n (T)dt, (58) 

0 

is sampled in t and T with a sampling period T, to obtain 

<7 A (0) = T[h( 0 + 0 )u n (0) + h( 0 + l)w„ (1) + • • • + h(0 + N)u n (N)] 

a „u n (1) = T[h{ 1 + 0)u„ (0) + h{\ + 1)M„ (1) + • • • + h(l + N)u n (N)] 


a n u n (N) = T[h(N + 0)u„ (0) + h(N + l)w„ (1) + • • • + h(N + N)u n ( W)] • (59) 

In this sampled approximation, the eigenvalues <7 n and eigenvectors u n (t ) are replaced by their 
approximate representations o n and u„(k) . Rewriting these equations in matrix form yields the 

eigenvalue relation 

a n [u n ]=TH[u n ], (60) 


where 



«„( 0 ) 

« n (D 


M N \ 


is the n 01 eigenvector and 6 n the corresponding eigenvalue of the matrix TH, with 

■/i(0) h( 1) ••• h(N) 

h( 1) h( 2) - h(N + l) 

H — 

h(N) h(N + l) ■■■ h(N + N) 


(61) 


(62) 


This eigenvalue problem may also be written as 

(< & n I-TH)[u n \ = 0 . (63) 

Instead of finding the Schmidt pairs as continuous functions of time, we may solve this 
eigenvalue problem to find them in sampled data form, as the eigenvectors of the sampled 
Hankel matrix H. The eigenvalues of H will approach the true Hankel singular values as the 
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sampling time gets small. Observe that the Hankel singular value approximations are the 
singular values of TH, not just H. The factor of T must be retained in order to maintain the 
correct magnitude of the resulting singular values. 

This approximation method has at least two sources of error. The most obvious error is 
in the sampling, or equivalently, the numerical evaluation of the integral. Although the error in 
the numerical approximation is readily obtainable, the relationship between this error and the 
resulting error in the Hankel singular values is not clear. Reducing the sampling time will 
increase the accuracy of the approximation, but it will also increase the size of the Hankel 
matrix, and will tend to make the matrix less well conditioned. The other obvious error is in the 
truncation of the integral. The limits of integration are from zero to infinity. Unfortunately, 
sampling in time out toward infinity will again yield more samples and a larger Hankel matrix. 
Thus this approximation method becomes a trade off. Decreasing T, or increasing N, improves 
the accuracy of the approximation but increases the computational complexity of the method. 
Again, the relationship between these time domain errors and the Hankel singular values and the 
Schmidt pairs is not clear. 

Returning to the norms related to the Hankel operator, it is now possible to approximate 
the norms by using the approximate Hankel singular values found from the above numerical 
technique. The magnitudes of the errors in this numerical process are not known. Nevertheless, 
these approximations are used in the next section to obtain system error bounds. 

Error Bounds based on the Approximate Hankel Operator 

It is now assumed that the signal of interest is the error signal obtained by differentiating 
the step response as in (27). This error signal represents the impulse response of the error 
system. As such, it can be used to create a Hankel matrix (62), after multiplication by the 
sampling time T. Consequently, all the results derived so far in this section may be applied to the 
modeling error. These results are repeated here in terms of t his error. 

The Hankel norm of the modeling error is defined to be the largest Hankel singular value, 
which is approximated by the largest singular value of the Hankel matrix of the sampled error 
impulse response; that is, 

||E(4,£<J, • (« 4 ) 
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The nuclear Hankel norm is then approximated by the sum of the singular values of the Hankel 
matrix of the error impulse response; that is, 

• < 65 > 

n=l 

The Hilbert-Schmidt Hankel norm is approximately the square root of the sum of the squares of 
the singular values of the sampled Hankel matrix; that is, 

( v V ' 2 

l £ < s >t- s -» = £ rf . 2 • (66) 

l n=1 

The Hilbert-Schmidt Hankel norm can also be approximated by the sampled time-weighted 
impulse energy using the relation 

< 67 > 

n=l Jk=0 

The following relationships can then be applied to the various error bounds 

\\ e ( S )\\ h < |£Og>)|L ^ HI, ^ 2 II £ ( s )L_ w • ( 6g ) 

The Hankel norm provides a lower bound on the frequency-domain infinity-norm of the error, 
while the time-domain 1-norm provides an upper bound. This means that the error system 
Nyquist plot is certain to lie entirely within a circle of radius ||e|| , but is also certain to protrude 
from a circle of radius ||E(j)|| for some frequencies. Twice the nuclear Hankel norm is an 
alternative upper bound on ||E(_/a)|| M ; so, comparing ||E(s)|| w , , and 2||£(s)||^_ w can provide 

some measure of the conservativeness of the error bounds. 

The Hilbert-Schmidt Hankel norm provides some independent insight into the quality of 
the model. Recall, for example, that |E(s)f _ ;/ is the area enclosed by the Nyquist plot of the 
error system. Thus, the ratio ||E(j)||^_ s _ w /lH(s)f H S H can be thought of as a measure of the 
relative Nyquist-plane modeling error. 
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6. MIMO Bounds 


Extension of the SISO Results 


The approach to the analysis of the modeling errors for multiple-input-multiple-output 
(MIMO) systems is similar to that described for single-input-single-output (SISO) systems in the 
previous sections. It is necessary first to define the MIMO error system by comparing the step- 
response matrix of a truth model with that of an approximate linear model. The difference 
between the step-response matrices is the step-response matrix of the error system. As before, 
this matrix of errors may be differentiated numerically to obtain the impulse-response matrix of 
the error system. The impulse-response matrix is written as 


e(t) = 


e \\ (0 

e 12 (o ••• 


^21 (0 

^22 (0 


.* rl (0 

e rl (t) ••• 

0. 


(69) 


where r is the number of outputs and m is the number of inputs. The Laplace transform of e(t) is 
the unknown error system transfer-function matrix E(s). 

In order to evaluate the accuracy of the MIMO model and hence its suitability for use in 
control system design, it is desirable to approximate the frequency-domain infinity-norm of 
E(jco). This is defined as 

\\EU(»)\l = maxb^ {£(;«)}] , (70) 

where <7^ {E(jco)} denotes the largest singular value of E(j(o). Thus the infinity-norm is the 
maximum value over all frequencies of the largest singular value of the error system matrix. The 
maximum singular value represents the largest possible ratio of the output vector magnitude to 
the input vector magnitude at each frequency. For the SISO case, the definition (70) reduces to 
the previously given definition (9), as the maximum singular value simply reduces to the 
magnitude of E(jco). 

A bound on ||E( i /co)|| oo may be obtained from the time-domain impulse-response error 
matrix according to the inequality 
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1M4 

IM»| • 

•• Ik (0||, 

|£(yn»L<CT nux 

IK, (oil, 

|^22 (0|j 

|Km(*)||, 


JKiCOlL 

IM4 • 

•• IM4 


The 1-norm of each element of the matrix on the right can be evaluated numerically as in the 
SISO case. The maximum singular value of this matrix is then the bound for the frequency- 
domain infinity-norm. The inequality (71) is found by extension of the relation (23) as follows: 
Let the notation Abs[A] denote the matrix whose elements are the magnitudes of the elements of 
the matrix A. Then, for any frequency 0), 

<*™x {E(jOi)}= <7,™ | j e(t) exp(-jox)dt^ 

OO OO 

£ CT^ | Abs[e(t ) exp(-jox j\dt = cr^ J Abs [e{t)]dt , (72) 

0 0 

where the final integral is exactly the same as the matrix on the right-hand side of (71). Since 
(72) holds for any (O, it also holds when cr^iEijay)} is at its maximum value \\E(j03)\\^ ; hence, 
(71) holds. 

MEMO system norms based on the Hankel operator are defined in the same way as in the 
SISO case. To obtain the approximate Hankel singular values, it is necessary to construct a 
block Hankel matrix containing the sampled error system impulse response. The error at each 
sampling instant will be an rxm matrix, which will constitute one block of the Hankel matrix. 
The singular values of this Hankel matrix can then be used to compute approximations to any of 
the Hankel norms. As in the SISO case, the relations 


£ ll £ o'"»IL 5 2 II £ mIL„ ( 73) 

are still valid. 

Scaling of Variables 

The actual value of any system norm will depend on the units chosen for the input and 
output variables. Therefore, it is necessary to consider the scaling of the input and output 
variables in order to evaluate whether a given system norm is “large” or “small.” The scaling is 
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especially important for a MIMO system, because there is a need to compare the relative sizes of 
the various transfer functions that constitute the system. If the variables are poorly scaled, the 
effect of a single input or output may dominate the value of the system norm, and the other 
inputs and outputs may have almost no effect on it. 

Figure 1 1 illustrates the scaling of the input and output variables. The actual system (or 
truth model) is depicted as the central block, which relates the input vector u to the output vector 
y. At a given operating point, the equilibrium values of u and y are denoted w e q and and the 
perturbations are denoted du and 8y, so that 


U = Weq + du. 

(74) 

y = y eq + dy. 

(75) 


Before the scaling, the equilibrium values u eq and yeq are effectively removed from consideration 
by use of the two summing junctions that represent (74-75). This leaves the perturbations du and 
5y, which are to be related by the linear model, but which may be badly scaled. A scaling factor 
is included with each input and output perturbation as 


du = M in du , 

(76) 

fy = M om<5y • 

(77) 


The scaling factors are included in the diagonal matrices M m and M out , and are chosen so that a 
full-scale perturbation du or dy corresponds to a value of ±1 for the scaled variable du or dy . A 
full-scale perturbation may be taken as the largest one tha: would ordinarily occur, in the best 
estimation of the modeler. In some cases the full-scale perturbation is assumed to be the nominal 
equilibrium value, so that the scaled-variable values represent the percentage (or fractional) 
changes from the nominal values. In this case, the diagonal elements of Mj n and M out are chosen 
as the elements of u eq and y^. 

Once the appropriate scaling factors are chosen, the value of the system norm should be 
computed using the input and output data represented as di and dy . It is assumed that output 
data are collected for separate step changes applied at all the inputs. Generally, the amplitudes of 
the applied step changes do not correspond to unit-amplitude inputs du ; therefore, in calculating 
the norm, the outputs must be scaled by the reciprocals of the amplitudes of the scaled step 
inputs that produced them. 
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7. Summary of Error Bounding Procedure 

This section contains a summary of the error bounding procedure presented in the 
foregoing. 

Step 1. Using a truth model or measured data, determine the equilibrium values of the system 
inputs and outputs at the desired operating point. Obtain a “truth” small-perturbation 
step response, and calculate Su and Sy, the input and output perturbations from the 
equilibrium. 

Step 2. Scale the small-perturbation variables Su and <5y by the respective chosen scaling 
factors A/" 1 and to obtain the non-dimensional perturbations Su and Sy . 

Step 3. Create a linear system model using some appropriate modeling method, and obtain its 
response Sy to the given scaled input perturbation Su . 

Step 4. Obtain the sampled step-error sequence e step ( kT ) = Sy{kT) - Sy(kT) . 

Step 5. Numerically differentiate the step response error to obtain the error-system impulse 
response. This can be done by means of a simple finite difference formula such as 

T 

Step 6. To obtain the normalized (unit) impulse response of the error system, divide e{kT) by 
the magnitude of Su . In the MEMO case, normalize each column of the impulse- 
response error matrix by the magnitude of the corresponding input Su i . The 
normalization could be performed instead on e slep (kT) in Step 4. 

Step 7. Approximate the 1-norm of the normalized error-system impulse response, 

\\e(t)l=jl T \e(k T )\. 

k = 0 

This is a bound for the frequency-domain error according to (26). It also gives a 
bound on the maximum output-error amplitude according to (29). 

a. The model frequency response can be plotted in the Nyquist plane along with 
the error-bound circle plotted at the origin, as in Figure 6. The model is 
relatively accurate at those frequencies for which the plot is outside this circle. 

b. The error bound can also be applied to the loop transfer function so as to 
analyze the stability robustness of a control design, as in Figures 8-10. 

In the MIMO case, compute the approximate norm of each element of the error 
matrix, and find the bound according to (71). The Nyquist-plane interpretation is 
difficult in this case. 

Step 8. The process can end here, or can continue using the Hankel operator. To pursue the 
Hankel-operator analysis, form the error-system Hankel matrix from the impulse- 
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response error. Find its singular values and left and right singular vectors as 
approximations to the Hankel singular values and Schmidt pairs, respectively. 

Step 9. From the Hankel singular value approximations in Step 8, find the Hankel norm, the 
nuclear norm, and the Hilbert-Schmidt norm of the error system. Circles with Hankel 
norm radius and twice the nuclear norm radius can also be drawn in the Nyquist plane 
to aid in understanding the conservativeness of the bounds. 

Step 10. If desired, determine the Hankel correction from the singular vector expansion in 
Step 8. Add the correction in parallel with the original system model to further 
reduce the modeling error and its bounds. 

This process assumes that the deviation of the linear small perturbation model from the 
truth model is roughly monotonic as measured by some norm of the error. This implies that the 
modeling error bounds are accurate as long as the system inputs are bounded by the original step 
input in their time-domain infinity norms. 

8. Example: A Mixed Compression Inlet 

This section contains an example of the calculation of error norms in the modeling 
process. The system chosen is a typical mixed-compression inlet (Chicatelli and Hartley, 1997), 
with the downstream corrected mass flow as the input, and the static pressure just downstream of 
the normal shock as the output. The truth-model data is taken from the nonlinear large- 
perturbation CFD simulation code LAPIN (Varner et al., 1987), which uses 123 nonlinear state- 
variable equations for this particular case. The linearization and reduction technique presented in 
(Chicatelli and Hartley, 1997) provided a lb^-order linear model. 

A constant mass-flow input u e q and a steady-state pressure output are taken from the 
LAPIN data. The scaling factor for the output perturbations is taken as 

M out = 117.8 lb/ft 2 , 

which is the free-stream pressure. The scaling factor for the input perturbations is taken as 

Mm = 1 Slug/S. 

Truth-model data were available for a l-% step change in the mass flow. The scaled 
(non-dimensional) step-response data are plotted in Figure 12 for both the truth model and the 
1 b^-order linear model. The difference between these two responses is the step response error. 
Dividing this error by the step-input magnitude (to normalize) and using (27) to estimate its 
derivative yields the error impulse response shown in Figure 13. The 1-norm of this signal is 
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|4 = 0.4633 , 

which is a bound on the frequency-domain infinity-norm of the error system. A circle of this 
radius, centered at the origin, is superimposed on the frequency-response polar plot of the 16 th - 
order linear model in Figure 14. This bound is also shown on the frequency-response magnitude 
plot as a straight line at -6.68 dB in Figure 15. Note that the system frequency response 
magnitude is larger than the error bound for all frequencies up to 100 Hz. This indicates that the 
model will be relatively accurate within that bandwidth. 

The error-system 1-norm also acts as a bound on the time-domain magnitude of any 
output error. According to (29), for small input variations (smaller than 1% in this case), the 
non-dimensional output error should remain less than 0.4633 times the non-dimensional input 
variation. 

The approximate Hankel singular values of the error system (i.e., the singular values of 
the Hankel matrix of the error sequence multiplied by the sampling interval T = 0.001 s) are 
plotted in Figure 16. According to (64), the largest of these singular values, 

=0.3249, 

represents the Hankel norm of the error system. This is a lower bound for the frequency-domain 
infinity-norm, meaning that the frequency-domain error magnitude is larger than this at some 
frequency. 

The nuclear Hankel norm is approximated using (65) as 

=4.0688. 

n = 1 

This relatively large number essentially means that any Hankel correction would have to be of 
high order to accomplish any significant further reduction of error. This can also be seen from 
the plot of the Hankel singular values in Figure 16, which reveals that cutting the Hankel-norm 
error in half would require at least a sixth-order correction. 

The Hilbert-Schmidt Hankel norm (or Frobenius Hankel norm) is approximated using 

(66) as 
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Recall that the area enclosed by the Nyquist plot of the error system is equal to the square of this 
Hilbert-Schmidt Hankel norm, or 0.7130. By comparison, the area enclosed by the Nyquist plot 
of the system model H(s ) itself is found to be 0.9018. 

9. Conclusion 

The step response of a system may be readily available, even when the impulse response 
and the frequency-domain description are not. This paper has shown how some time-domain 
and frequency-domain bounds useful in modeling and control design may be determined from 
the system step response. The calculation of the bounds requires first a numerical differentiation 
of the step-response error to estimate the impulse-response error. A frequency-domain error 
bound may then be determined by use of a simple numerical integration. The Hankel norms of 
the error system may also be determined from a Hankel matrix of the impulse-response error, by 
use of singular-value calculations. The Hankel-norm information may be used to compute some 
additional frequency-domain upper and lower bounds, to determine an error-reducing correction 
to a given model, and to qualitatively evaluate the relative frequency-domain modeling error in 
the Nyquist plane. 
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Appendix: Matlab code listing for computing example error bounds 

% FILE norms . m 
% 

% Data file 'normdata' contains scaled output data from the truth 
% model (ynl), and unsealed output data from the 16th- order model (y 

load normdata 
T=t ime ( 10 ) -time ( 9 ) ; 

% 

% Plot truth-model and 16th~order model step responses 

figure ( 12 ) 

plot (time , ynl ( : , 1 ) , 1 o ' , t ime , ynl (1,1) +yr ( : , 1 ) ) 
xlabel ( ' Time (seconds)’) 

ylabel ( 1 Scaled (non-dimensional ) pressure ' ) 

text (0.953,10.03, ’ o truth model ' ) 

text ( 0 . 953 , 10 . 01 , ’ - 16th-oraer linear model') 

grid 

CL 

t* 

% Calculate scaled (non-dimensional ) error , 

e= (ynl ( : , 1) -ynl (5,1) *ones (size(yr ( :,!)))) -yr( : , 1) ; 

% 

% Differentiate to find impulse- response error. 

ei=diff (e) /T; 

% Divide by step- input magnitude to normalize ei 

ei=ei/ ( (input (190 ) -input (5) ) ) ; 

% 

% Plot scaled, normalized inpulse- response error. 

figure ( 13 ) 

plot (time (1:196), ei, * b ' ) 

ylabel ( ' Scaled, normalized impulse-response error 1 ) 

xlabel ( 'Time (seconds) ' ) 

grid 

a 

G 

% Compute one -norm of the impulse- response error. 

onenorm=sum ( abs ( ei ) ) *T 

9 - 

% Calculate coordinates of circle in complex plane. 

j=sqrt (-1) ; 
w= [-180:180] *pi / 180 ; 
cir=onenorm*exp ( j *w) ; 

% 

% Cal cu 1 ate f r e qu eric y res pons e of 1 6 1 h - o r cle r mode 1 . 

sysr=ss (ar , br , cr ( 1 , : ) , dr ( 1 , : ) ) ; 
wr=logspace (0,4,200) ; 

[rfr , ifr] =nyquist (sysr, wr ) ; 

% Scale the frequency response: 117,8 is the freestream pressure 
% This scaling was not built into the 16th- order model 

rf r=rf r/117 . 8 ; 
if r=if r/117 . 8 ; 

% 

% Make polar plot with circle to show error bound. 

figure ( 14 ) 
plot (rfr , if r ) 
axis equal 
hold on 

plot (real (cir) , imag(cir) , 'r ’ ) 
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ylabel (' Imaginary H(jw)*) 
xlabel ('Real H(jw) 1 ) 
grid 
hold off 

% 

% Make magnitude plot with line to show error bound. 

figure ( 15 ) 

wr=logspace (0,3,200) ; 

[rfr, ifr] =nyquist ( sysr, wr ) ; 
rf r=rf r/117 . 8; 
ifr=ifr/117 .8; 

semilogx (wr/6 .28, 20*1 og 10 ( sqrt (rfr.*rfr+ifr.*ifr) ) ) 
hold on 

semilogx (wr/6 .28,20*logl0 ( onenorm*ones (size(rfr) ) ) , ' r 1 ) 
xlabel (' Frequency (Hz)') 

ylabel ( 1 Frequency -response magnitude (dB) ' ) 
grid 

hold off 

% 

% Construct Ilankel matrix from impulse- response error sequence. 

h=hankel ( ei ) ; 

% Truncate matrix to eliminate zeros. 

m= length ( ei ) / 2 ; 
h=h ( 1 : m, 1 : m) ; 

% Calculate Hankel singular values 

[u, s , v] =svd ( h ) ; 

% Put Hankel singular values into a vector 

s=diag ( s*T) ; 

% 

% Plot Ilankel singular values. 

figure ( 16 ) 
plot ( s , ’o') 
grid 

ylabel (' Hankel singular value') 
xlabel (' Singular value number') 

% 

% Calculate Hankel norms. 

hnorm=s ( 1 ) 
nnorm=sum ( s ) 
hsnorm=sqrt ( sum ( s . *s ) ) 
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Figure 3. Frequency-domain error bound represented in the Nyquist plane. 
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Figure 4. Obtaining the system impulse response the hard way. 
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Figure 7. Closed-loop system with unmodeled dynamics included. 
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Figure 8. Frequency-dependent error bounds represented in the Nyquist plane. 
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Figure 9. Worst-case error bound represented in the Nyquist plane. 
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Figure 10. Worst-case error bound represented in the Nyquist plane, revisited. 
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Figure 11. Scaling of the input and output variables. 
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Figure 13. Impulse response error found by numerical differentiation of the step response error. 



Figure 14. Polar plot of frequency response of 1 b^-order model. 
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Figure 15. Magnitude plot of frequency response of 1 6 lh -order model. 



Figure 16. Hankel singular values of the error system. 
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